run scripts/start_Smithkline.m

%% Import data
opts = spreadsheetImportOptions("NumVariables", 6);
addpath('figures')

% Specify sheet and range
opts.Sheet = "ForMatlab_Final";
opts.DataRange = "A1:F1699";

% Specify column names and types
opts.VariableNames = ["DateVec", "BSKHPrice", "BSKH2Price", "BSKHVolume", "BSKEPrice", "BSKEVolume"];
opts.VariableTypes = ["string", "double", "double", "double", "double", "double"];

% Specify variable properties
opts = setvaropts(opts, "DateVec", "WhitespaceRule", "preserve");
opts = setvaropts(opts, "DateVec", "EmptyFieldRule", "auto");

% Import the data
tbl = readtable("raw_data\SmithKline data.xlsx", opts, "UseExcel", false);

%% Convert to output type
DateVec    = tbl.DateVec;
BSKHPrice  = tbl.BSKHPrice;
BSKH2Price = tbl.BSKH2Price;
BSKHVolume = tbl.BSKHVolume;
BSKEPrice  = tbl.BSKEPrice;
BSKEVolume = tbl.BSKEVolume;

clear opts tbl

DateVec    = datetime(DateVec);
Dates      = DateVec;
BSKHPrice  = BSKHPrice(1:end,:);
BSKHVolume = BSKHVolume(1:end,:);
BSKEPrice  = BSKEPrice(1:end,:);
BSKEVolume = BSKEVolume(1:end,:);

meanDVs = zeros(size(Dates,1),2);
counter = 0;
while (counter<size(Dates,1))
   counter = counter + 1; 
   if(counter<261)
       fracDVs(counter,:) = [0 0];
   else
       %Multiply by 10 because of scaling: prices are rescaled up by 100
       %and volumes are scaled down by 1000 so 1000/100 = 10 for true value
       fracDVs(counter,:) = [mean(BSKHVolume(counter-260:counter,1).*(BSKHPrice(counter-260:counter,1))) mean(BSKEVolume(counter-260:counter,1).*(BSKEPrice(counter-260:counter,1)))];
   end
end

arbSizes = zeros(size(Dates,1)-260,1);
totalWelfare = zeros(size(Dates,1)-260,1);
fnEvals = zeros(size(Dates,1)-260,1);

counter = 260;
while(counter<size(Dates,1))
   counter = counter + 1;
   Values = @(m) (2*BSKHPrice(counter,1).*(1+(100*(1000000*m/fracDVs(counter,1))*0.0000) + (((100*abs(1000000*m)./fracDVs(counter,1)).^0.5)*0.000889.*sign(m))))-(BSKEPrice(counter,1)*(1-(100*(1000000*m/fracDVs(counter,2))*0.0000) - (((100*abs(1000000*m)/fracDVs(counter,2)).^0.5)*0.000889.*sign(m))));
   Values2 = @(m) (2*BSKHPrice(counter,1).*(1+(100*(1000000*m/fracDVs(counter,1))*0.0000) + (((100*abs(1000000*m)./fracDVs(counter,1)).^0.5)*0.000889.*sign(m))))./(BSKEPrice(counter,1)*(1-(100*(1000000*m/fracDVs(counter,2))*0.0000) - (((100*abs(1000000*m)/fracDVs(counter,2)).^0.5)*0.000889.*sign(m))))-1;
   
   options = optimoptions('fsolve','FunctionTolerance',1e-12,'OptimalityTolerance',1e-12,'StepTolerance',1e-6,'MaxFunctionEvaluations',10000);
   [m1,fnVal] = fsolve(Values,1,options);
   grid = sign(m1)*(0:1:abs(m1));
   evals2 = zeros(size(grid,2),1);
   counter2 = 0;
   while(counter2<size(grid,2))
       counter2 = counter2+1;
       evals2(counter2,1) = Values2(grid(1,counter2));
   end
   arbSizes(counter-260,1) = m1*1000000;
   %evals = Values(grid);
   totalWelfare(counter-260,1) = sum(evals2)*1000000;
   fnEvals(counter-260,1)=fnVal;
end

ratios = 2*BSKHPrice./(BSKEPrice);
ratios = ratios(261:end,:);

save('intermediate\workspace_Smithkline');

%% Create figures
load('intermediate\workspace_Smithkline');

fig_BSK_welfare(DateVec, ratios, totalWelfare, optfig)
fig_BSK_arbsize(DateVec, ratios, arbSizes, optfig)